################################################################################
## "Do Conspiracy Theories Undermine Support for Democracy? 
## Experimental Evidence from Brazil" (Amaral et al.)
## Online Support Information (OSI):
################################################################################

# This code:
# 1. Tests exploratory hypotheses
# 2. Generates tables and graphs in the Supporting Information

library(tidyverse)
library(haven)
library(sjlabelled)
library(psych)
library(sjPlot)
library(Hmisc)
library(lavaan)
library(broom)
library(stargazer)
library(Rmisc)

################################################################################
# C Pilot --------------------------------------------------------------------
################################################################################

# Data --------------------------------------------------------------------

setwd("C:/Users/vanes/Downloads")

## Pilot 1 dataset
p1 <- read_sav("C:/Users/vanes/Downloads/piloto.sav")
p1_csv <- read.csv("C:/Users/vanes/Downloads/piloto.csv")

## Pilot 2 dataset
p2 <- read_sav("C:/Users/vanes/Downloads/piloto_2.sav")

# Data cleaning -----------------------------------------------------------

## Changes 99 to NA in the pilot 1 data
p1[p1 == 99] <- NA

# Renames variables
p1 <- p1 %>% 
  mutate(affect_pt = P2_1,
         affect_bs = P2_2, 
         consp_left = P3,
         consp_right = P4,
         consp_freem = P5,
         belief = P10,
         consp_stab = P9_1,
         consp_gen = P9_2,
         polgroup = P8_1,
         ideolgroup = P8_2,
         exp_group = p1_csv$ROTACION,
         villain = p1_csv$P6,
         exposure = P11,
         victim = P7)

p2 <- p2 %>%
  mutate(affect_pt = P2_1,
         affect_bs = P2_2, 
         belief = P10,
         consp_stab = P9_1,
         consp_gen = P9_2,
         exposure = P11)

## Conspiracism index
p2 <- p2 %>%
  mutate(consp_index = (P9_1 + P9_2) / 2)

## Renames experimental groups
p2 <- p2 %>%
  mutate(exp_group = factor(ifelse(ROTACION == 1, 'Placebo', 'Control')))

table(p2$exp_group)

## Recodes categories of political affect variables, reverses conspiracionism
## and exposure measures, relevels experimental groups
p1 <- p1 %>%
  mutate(affect_bs_cat = factor(ifelse(affect_bs <= 3, 'Contempt',
                                       ifelse(affect_bs >= 4 & affect_bs <= 6, 'Neutrality',
                                              ifelse(affect_bs >= 7, 'Affection', NA))), 
                                levels = c('Contempt', 'Neutrality', 'Affection')),
         affect_pt_cat = factor(ifelse(affect_pt <= 3, 'Contempt',
                                       ifelse(affect_pt >= 4 & affect_pt <= 6, 'Neutrality',
                                              ifelse(affect_pt >= 7, 'Affection', NA))), 
                                levels = c('Contempt', 'Neutrality', 'Affection')),
         belief = max(belief, na.rm = T) - belief,
         consp_left = max(consp_left, na.rm = T) - consp_left,
         consp_right = max(consp_right, na.rm = T) - consp_right,
         consp_freem = max(consp_freem, na.rm = T) - consp_freem, 
         consp_gen = max(consp_gen, na.rm = T) - consp_gen,
         consp_stab = max(consp_stab, na.rm = T) - consp_stab,
         exp_group = factor(exp_group,
                            levels = c('Teoria Oficial', 'Teoria Maconaria',
                                       'Teoria Left Behind', 'Teoria Fake Stabbing')),
         exposure = max(exposure, na.rm = T) - exposure,
         polgroup = max(polgroup, na.rm = T) - polgroup)



# Table A1 ----------------------------------------------------------------

## Regression models predicting belief 
summary(
  lm(
    belief ~ affect_pt + affect_bs,
    data = p1 %>% filter(exp_group == 'Teoria Maconaria')) -> belief1
)

summary(
  lm(
    belief ~ affect_pt + affect_bs,
    data = p1 %>% filter(exp_group == 'Teoria Fake Stabbing')) -> belief2
)

summary(
  lm(
    belief ~ affect_pt + affect_bs,
    data = p1 %>% filter(exp_group == 'Teoria Left Behind')) -> belief3
)

summary(
  lm(
    belief ~ affect_pt + affect_bs,
    data = p1 %>% filter(exp_group == 'Teoria Oficial')) -> belief4
)

summary(
  lm(
    belief ~ affect_pt + affect_bs,
    data = p2 %>% filter(exp_group == 'Placebo')) -> belief5
)

## Generates table with stargazer
stargazer(belief5, belief4, belief1, belief3, belief2, 
          column.labels = c('Placebo',
                            'Official',
                            'Freemasonry',
                            'Left',
                            'Right'), 
          covariate.labels = c('Affect Petistas', 
                               'Affect Bolsonaristas'), 
          dep.var.labels = 'Belief')


# Table A2 ----------------------------------------------------------------

## OLS models predicting belief in different groups as responsible for Bolsonaro's
## attempted murder

summary(
  lm(
    consp_freem ~ affect_pt_cat,
    data = p1) -> freembehind
)

# is the pt behind?
summary(
  lm(
    consp_left ~ affect_pt_cat,
    data = p1) -> ptbehind
)

# is bolsonaro behind?
summary(
  lm(
    consp_right ~ affect_pt_cat,
    data = p1) -> rightbehind
)

# table with coefficients
stargazer::stargazer(freembehind, ptbehind, rightbehind,
                     dep.var.labels = c('Freemasonry', 'Political Left and PT', 
                                        'Political Right and Bolsonaro'),
                     covariate.labels = c('PT Neutrality', 'PT Affection'))


# Table A3 ----------------------------------------------------------------

## Dumies for perpetrator groups
p1 <- p1 %>%
  mutate(left_pt = ifelse(villain == 'A esquerda e o PT', 1, 0),
         right_bl = ifelse(villain == 'A direita e Bolsonaro', 1, 0),
         freem = ifelse(villain == 'A maçonaria', 1, 0),
         mourao = ifelse(villain == 'Hamilton Mourão', 1, 0),
         bispo = ifelse(villain == 'Adélio Bispo', 1, 0),
         novillain = ifelse(villain == 'Não é possível identificar um vilão', 1, 0))

## OLS regressions predicting the group perceived as the perpetrator

vil1 <- lm(data = p1, 
           left_pt ~ exp_group)

vil2 <- lm(data = p1, 
           right_bl ~ exp_group)

vil3 <- lm(data = p1, 
           freem ~ exp_group)

vil4 <- lm(data = p1, 
           mourao ~ exp_group)

vil5 <- lm(data = p1, 
           bispo ~ exp_group)

vil6 <- lm(data = p1, 
           novillain ~ exp_group)

## Generates table A3
stargazer(vil1, vil2, vil3, vil4, vil5, vil6, type = 'latex',
          covariate.labels = c('Freemasons',
                               'Left Behind',
                               'Fake Stabbing'),
          dep.var.labels = c('Left/PT', 
                             'Right/Bolsonaro', 
                             'Freemasons', 
                             'Mourão', 
                             'Bispo',
                             'No Villain'), 
          omit.stat = c('f', 'ser'))


# Table A4 ----------------------------------------------------------------

## Dummies for victimized group
p1 <- p1 %>%
  mutate(v_left_pt = ifelse(victim == 3, 1, 0),
         v_bols = ifelse(victim == 2, 1, 0),
         v_right = ifelse(victim == 1, 1, 0),
         v_freem = ifelse(victim == 4, 1, 0),
         v_mourao = ifelse(victim == 5, 1, 0),
         v_bispo = ifelse(victim == 6, 1, 0),
         v_novictim = ifelse(victim == 7, 1, 0))

vic1 <- lm(data = p1, 
           v_left_pt ~ exp_group)

vic2 <- lm(data = p1, 
           v_right ~ exp_group)

vic3 <- lm(data = p1, 
           v_freem ~ exp_group)

vic4 <- lm(data = p1, 
           v_mourao ~ exp_group)

vic5 <- lm(data = p1, 
           v_bispo ~ exp_group)

vic6 <- lm(data = p1, 
           v_novictim ~ exp_group)

vic7 <- lm(data = p1, 
           v_bols ~ exp_group)

stargazer::stargazer(vic1, vic7, vic2, vic3, vic4, vic5, vic6,
                     covariate.labels = c('Freemasons',
                                          'Left Behind',
                                          'Fake Stabbing'),
                     dep.var.labels = c('Left/PT', 'Bolsonaro', 'Right', 
                                        'Freemasons', 'Mourão', 'Bispo',
                                        'No Victim'),
                     omit.stat = c('f', 'ser'))


# Figure A1 ---------------------------------------------------------------

## Generates graph
consp_groups_level <- p1 %>%
  select(consp_left, consp_right, consp_freem) %>%
  pivot_longer(cols = c(consp_left, consp_right, consp_freem), 
               names_to = 'Variable', values_to = 'val') %>%
  summarySE(measurevar = 'val', groupvars = 'Variable', na.rm = T) %>%
  ggplot(aes(x = Variable, y = val)) + 
  geom_errorbar(aes(ymin = val - ci, ymax = val + ci), width = .2) +
  geom_point() +
  theme_bw() +
  labs(y = 'Potential conspirator') +
  theme(plot.caption = element_text(hjust = 0),
        plot.title = element_text(hjust = .5), 
        panel.grid = element_blank())

consp_groups_level

## Exports figure 
ggsave('./potentialconspirators.png', plot = consp_groups_level, 
       device = 'png', width = 3, height = 4)



# Figure A2 ---------------------------------------------------------------


## Creates conspiracism index
p1 <- p1 %>%
  mutate(consp_index = (consp_gen + consp_stab) / 2)

p2 <- p2 %>%
  mutate(consp_index = (P9_1 + P9_2) / 2)

## Experimental group labels
#p1 <- p1 %>%
#mutate(exp_group = recode_factor(df$exp_group, 
#`Teoria Oficial` = 'Official Theory',
#`Teoria Maconaria` = 'Freemasons Theory',
#`Teoria Left Behind` = 'Left Behind Theory',
#`Teoria Fake Stabbing` = 'Fake Stabbing Theory'))

p1 <- p1 %>%
  mutate(
    exp_group = recode_factor(
      exp_group,
      `Teoria Oficial` = "Official Theory",
      `Teoria Maconaria` = "Freemasons Theory",
      `Teoria Left Behind` = "Left Behind Theory",
      `Teoria Fake Stabbing` = "Fake Stabbing Theory"
    )
  )

## Generates graph
conspiracism <- p1 %>%
  select(consp_index, exp_group) %>%
  mutate(Pilot = 'Pilot 1') %>%
  bind_rows(p2 %>% 
              select(consp_index, exp_group) %>%
              mutate(Pilot = 'Pilot 2')) %>%
  mutate(exp_group = factor(exp_group, 
                            levels = c('Control', 
                                       'Placebo',
                                       'Official Theory', 
                                       'Freemasons Theory',
                                       'Left Behind Theory',
                                       'Fake Stabbing Theory'))) %>%
  mutate(exp_group = recode_factor(exp_group, 
                                   Control = 'Control',
                                   Placebo = 'Placebo',
                                   `Official Theory` = 'Official',
                                   `Freemasons Theory` = 'Freemasons',
                                   `Left Behind Theory` = 'Left Behind',
                                   `Fake Stabbing Theory` = 'Fake Stabbing')) %>%
  summarySE(measurevar = 'consp_index', 
            groupvars = c('exp_group', 'Pilot'), na.rm = T) %>%
  ggplot(aes(x = exp_group, y = consp_index, color = Pilot)) + 
  geom_errorbar(aes(ymin = consp_index - ci, ymax = consp_index + ci), 
                width = .2, position = position_dodge(width = 0.5)) +
  geom_point(position = position_dodge(width = 0.5)) +
  theme_bw() +
  labs(y = 'Levels of Conspiracism', x = 'Experimental Group') +
  theme(plot.caption = element_text(hjust = 0),
        plot.title = element_text(hjust = .5), 
        panel.grid = element_blank(),
        axis.text.x = element_text(angle = 40, hjust = 1, vjust = 1))

conspiracism

ggsave('./conspiracism.png', conspiracism, 'png',
       width = 5, height = 4)



# Figure A3 ---------------------------------------------------------------


## Generates graph

belief_exposure <- p1 %>%
  select(exposure, belief, exp_group) %>%
  mutate(Pilot = 'Pilot 1') %>%
  bind_rows(p2 %>% 
              select(exposure, belief, exp_group) %>%
              mutate(Pilot = 'Pilot 2')) %>%
  mutate(`Previous Exposure` = exposure,
         `Belief` = belief) %>%
  select(`Previous Exposure`, Belief, exp_group, Pilot) %>%
  pivot_longer(cols = !c(exp_group, Pilot), names_to = 'variable', values_to = 'value') %>%
  filter(exp_group != 'Control') %>%
  mutate(exp_group = factor(exp_group, 
                            levels = c('Placebo',
                                       'Official Theory', 
                                       'Freemasons Theory',
                                       'Left Behind Theory',
                                       'Fake Stabbing Theory'))) %>%
  mutate(exp_group = recode_factor(exp_group, 
                                   Placebo = 'Placebo',
                                   `Official Theory` = 'Official',
                                   `Freemasons Theory` = 'Freemasons',
                                   `Left Behind Theory` = 'Left',
                                   `Fake Stabbing Theory` = 'Right')) %>%
  summarySE(measurevar = 'value', 
            groupvars = c('exp_group', 'variable', 'Pilot'), na.rm = T) %>%
  ggplot(aes(x = exp_group, y = value, color = Pilot)) + 
  geom_errorbar(aes(ymin = value - ci, ymax = value + ci), 
                width = .2, position = position_dodge(width = 0.5)) +
  geom_point(position = position_dodge(width = 0.5)) +
  theme_bw() +
  labs(y = '', x = 'Experimental Group') +
  theme(plot.caption = element_text(hjust = 0),
        plot.title = element_text(hjust = .5), 
        panel.grid = element_blank(),
        axis.text.x = element_text(angle = 40, hjust = 1, vjust = 1)) +
  facet_wrap(~variable)

belief_exposure

ggsave('./belief_exposure.png', belief_exposure, 'png',
       width = 6, height = 4)

################################################################################
# D Study 1 - Supplementary Analyses ------------------------------------------
################################################################################

rm(list = ls())

# D.1 Balance Tests --------------------------------------------------

## Study 1 dataset
df1 <- read.csv('./Study 1_clean_pca outcomes.csv')

## Study 2 dataset
df2 <- read.csv('./Study 2_clean_pca outcomes.csv')

## Changes reference level of the experimental group variable
df1$exp_group <- factor(df1$exp_group, levels = c('Placebo',
                                                  'Oficial',
                                                  'Maçonaria',
                                                  'Esquerda',
                                                  'Direita'))

df2$exp_group <- factor(df2$exp_arm, levels = c('Placebo',
                                                'Official',
                                                'Freemasonry',
                                                'Left',
                                                'Right'))

# Table A5 ----------------------------------------------------------------

## Changes variables to numeric format
df1 <- df1 %>%
  dplyr::mutate(pol_interest = as.numeric(pol_interest),
                affect_bs = as.numeric(affect_bs),
                affect_lula = as.numeric(affect_lula),
                lr = as.numeric(lr),
                education = as.numeric(F5),
                sex = as.numeric(F1),
                income = as.numeric(F9),
                age2 = as.numeric(F2))


df1 <- df1 %>%
  mutate(exp_group_en = recode_factor(
    exp_group,
    Placebo   = "Placebo",
    Oficial   = "Official",
    Maçonaria = "Freemasonry",
    Esquerda  = "Left",
    Direita   = "Right"
  ))

vtable::st(
  df1,
  vars = c("pol_interest","affect_bs","affect_lula","lr","education","sex","age2","income"),
  group = "exp_group_en",
  group.test = TRUE,
  out = "latex",
  file = "table_A5_sumstats_study1.tex",
  digits = 4,
  title = "Balance Tests"
)


# testing relationship between pre-treatment variables and compliance
summary(
  lm(scale(compliance_n) ~ pol_interest + affect_bs + affect_lula + lr + education + 
       sex + age + income,
     df1)
)

# D.2 Descriptive Statistics of Outcomes ----------------------------------------

# Table A6 ----------------------------------------------------------------

## Creates data.frame with outcome variables
outcomes <- df1[, c('acc_results',
                    'acc_results_cand',
                    'strongman',
                    'pol_influence',
                    'ban_protests',
                    'ban_critics',
                    'violence',
                    'sup_dem',
                    'agression_bs',
                    'agression_pt')] 

## Changes column names
names(outcomes) <- c('Accept election results',
                     'Candidates should accept election results',
                     'Support for strongmen',
                     'Undue political influence',
                     'Ban protests',
                     'Ban government critics',
                     'Electoral violence',
                     'Support for democracy',
                     'Aggression bolsonaristas',
                     'Aggresion lulistas')

## Generates table
stargazer(as.data.frame(outcomes), type = 'latex', mean.sd = TRUE, 
          nobs = FALSE, median = FALSE, iqr = FALSE,
          digits=1, align=T, title = 'Summary Statistics')

# D.3 Robustness Checks ----------------------------------------

# Table A7 ----------------------------------------------------------------

df1 <- df1 %>%
  dplyr::mutate(complier = ifelse(exp_group == 'Placebo' & TIME_VIDEO > 110 |
                                    exp_group == 'Oficial' & TIME_VIDEO > 73 |
                                    exp_group == 'Maçonaria' & TIME_VIDEO > 124 |
                                    exp_group == 'Esquerda' & TIME_VIDEO > 104 |
                                    exp_group == 'Direita' & TIME_VIDEO > 118, 1, 0),
                attentive = ifelse(exp_group == 'Placebo' & TIME_VIDEO > 2*110 |
                                     exp_group == 'Oficial' & TIME_VIDEO > 2*73 |
                                     exp_group == 'Maçonaria' & TIME_VIDEO > 2*124 |
                                     exp_group == 'Esquerda' & TIME_VIDEO > 2*104 |
                                     exp_group == 'Direita' & TIME_VIDEO > 2*118, 0, 1))

# splitting time_video variable into bins of 3 seconds each
df1$time_video_bin <- as.numeric(cut2(df1$TIME_VIDEO,
                                      cuts = seq(min(df1$TIME_VIDEO),
                                                 max(df1$TIME_VIDEO), 3)))

# creating function to calculate statistical mode
getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

## creating continuous variable for compliance
## here we create a weighting variable based on the inverse of the absolute
## difference from the modal seconds bin to the actual seconds bin of each
## observation

# we sum the absolute difference by one not to divide by zero
df1 <- df1 %>%
  group_by(exp_group) %>%
  dplyr::mutate(compliance_n = 1 / (abs(getmode(time_video_bin) - time_video_bin) + 1))


## Creates indicator to measure whether people drop out during the video
df1 <- df1 %>%
  group_by(exp_group) %>%
  mutate(attentive2 = ifelse(TIME_VIDEO <= quantile(TIME_VIDEO, .75) + 1.5*IQR(TIME_VIDEO), 0, 1))

summary(
  compliance <- lm(complier ~ exp_group,
                   df1)
)

stargazer(compliance, type = 'latex',
          dep.var.labels = 'Attrition',
          covariate.labels = c('Official', 'Freemasonry', 'Left', 'Right'),
          omit = c('Constant'),
          omit.stat = c('adj.rsq', 'f', 'ser'),
          add.lines = list(c('Baseline', 'Placebo')))



# Table A8 ----------------------------------------------------------------

## Study 1 dataset
df1 <- read.csv('./Study 1_clean_pca outcomes.csv')

df1 <- df1 %>%
  dplyr::mutate(complier = ifelse(exp_group == 'Placebo' & TIME_VIDEO > 112.75 |
                                    exp_group == 'Oficial' & TIME_VIDEO > 76.67 |
                                    exp_group == 'Maçonaria' & TIME_VIDEO > 126.9 |
                                    exp_group == 'Esquerda' & TIME_VIDEO > 107.3 |
                                    exp_group == 'Direita' & TIME_VIDEO > 121.25, 1, 0),
                attentive = ifelse(exp_group == 'Placebo' & TIME_VIDEO > 2*110 |
                                     exp_group == 'Oficial' & TIME_VIDEO > 2*73 |
                                     exp_group == 'Maçonaria' & TIME_VIDEO > 2*124 |
                                     exp_group == 'Esquerda' & TIME_VIDEO > 2*104 |
                                     exp_group == 'Direita' & TIME_VIDEO > 2*118, 0, 1))

table(df1$exp_group[df1$complier == 1])

# splitting time_video variable into bins of 3 seconds each
df1$time_video_bin <- as.numeric(cut2(df1$TIME_VIDEO, 
                                      cuts = seq(min(df1$TIME_VIDEO), 
                                                 max(df1$TIME_VIDEO), 3)))

# creating function to calculate statistical mode
getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

# creating continuous variable for compliance
# here we create a weighting variable based on the inverse of the absolute
# difference from the modal seconds bin to the actual seconds bin of each
# observation

# we sum the absolute difference by one not to divide by zero
df1 <- df1 %>%
  group_by(exp_group) %>%
  dplyr::mutate(compliance_n = 1 / (abs(getmode(time_video_bin) - time_video_bin) + 1))


summary(
  h1.p.late <- lm(scale(dem_index) ~ CTxPlacebo, df1 %>% filter(complier == 1))
)

summary(
  h1.p.wls <- lm(scale(dem_index) ~ CTxPlacebo, df1 %>% filter(complier == 1), 
                 weights = compliance_n)
)

summary(
  h1.o.late <- lm(scale(dem_index) ~ XCTxOficial, df1 %>% filter(complier == 1))
)

summary(
  h1.o.wls <- lm(scale(dem_index) ~ XCTxOficial, df1 %>% filter(complier == 1), 
                 weights = compliance_n)
)

# Baseline Placebo
d_placebo <- df1 %>%
  dplyr::filter(exp_group %in% c("Placebo","Maçonaria","Esquerda","Direita")) %>%
  dplyr::mutate(CT = ifelse(exp_group == "Placebo", 0, 1))

# Baseline Official
d_official <- df1 %>%
  dplyr::filter(exp_group %in% c("Oficial","Maçonaria","Esquerda","Direita")) %>%
  dplyr::mutate(CT = ifelse(exp_group == "Oficial", 0, 1))

## ITT
h1.p.itt <- lm(scale(dem_index) ~ CT, data = d_placebo)
h1.o.itt <- lm(scale(dem_index) ~ CT, data = d_official)

## LATE (compliers)
h1.p.late <- lm(scale(dem_index) ~ CT,
                data = d_placebo %>% dplyr::filter(complier == 1))

h1.o.late <- lm(scale(dem_index) ~ CT,
                data = d_official %>% dplyr::filter(complier == 1))

## WLS
h1.p.wls <- lm(scale(dem_index) ~ CT,
               data = d_placebo %>% dplyr::filter(complier == 1),
               weights = compliance_n)

h1.o.wls <- lm(scale(dem_index) ~ CT,
               data = d_official %>% dplyr::filter(complier == 1),
               weights = compliance_n)

stargazer::stargazer(
  h1.p.itt, h1.p.late, h1.p.wls,
  h1.o.itt, h1.o.late, h1.o.wls,
  type = "latex",
  dep.var.labels = "Overall Support for Democracy",
  covariate.labels = "CT",
  column.labels = c("ITT","LATE","WLS","ITT","LATE","WLS"),
  omit.stat = c("adj.rsq", "f", "ser"),
  add.lines = list(
    c("Baseline","Placebo","Placebo","Placebo","Official","Official","Official")
  ),
  digits = 3
)

################################################################################
## E Study 2 - Supplementary Analyses  ----------------------------------------
################################################################################

## E.1 Balance Tests  ----------------------------------------------------------

# Table A9 ----------------------------------------------------------------

rm(list = ls())

## Study 2 dataset
df2 <- read.csv('./Study 2_clean_pca outcomes.csv')

## Changes reference level of the experimental group variable

df2$exp_group <- factor(df2$exp_arm, levels = c('Placebo',
                                                'Official',
                                                'Freemasonry',
                                                'Left',
                                                'Right'))
# Changes variables to numeric format
df2 <- df2 %>%
  dplyr::mutate(Gender = as.numeric(as.factor(Q.gender)),
                Education = as.numeric(as.factor(Q.education)),
                Income = as.numeric(as.factor(Q.renda)),
                `Conspiracionism 1` = as.numeric(Q.conspiracionismo1),
                `Conspiracionism 2` = as.numeric(Q.conspiracionismo2),
                `Political Interest` = as.numeric(as.factor(Q.interesse)),
                `Left-Right` = as.numeric(Q.ideology),
                `Vote` = as.numeric(as.factor(Q.previousvote)),
                `Manipulation` = exp_arm)

# Generates table
vtable::sumtable(data = df2, out = 'latex',
                 vars = c('Gender', 'Education', 'Income', 'Conspiracionism 1',
                          'Conspiracionism 2', 'Political Interest',
                          'Left-Right', 'Vote'),
                 group = 'Manipulation', group.test = T, title = '')

## E.2 Descriptive Statistics of Outcomes ------------------------------------

# Table A10  ----------------------------------------------------------------

outcomes2 <- df2 %>%
  select(Q.trust.parties : Q.trust.media, Q.trust_eleições, Q.invasão,
         extr.outgroup, dehum.outgroup, Q.identidade_imp, Q.aceitação1, Q.aceitação2)

names(outcomes2) <- c('Trust parties',
                      'Trust Congress',
                      'Trust Judiciary',
                      'Trust police',
                      'Trust media',
                      'Trust elections',
                      'Justify Jan 8 invasion',
                      'Extremism outgroups',
                      'Dehumanization outgroups',
                      'Importance of political identity',
                      'Candidates should accept election results',
                      'Accept election results')

stargazer(as.data.frame(outcomes2), type = 'latex', mean.sd = TRUE, 
          nobs = FALSE, median = FALSE, iqr = FALSE,
          digits=1, align=T, title = 'Outcome Descriptive Statistics (Study 2)')


# E.3 Main Models --------------------------------------------------------------

## Table A11 -------------------------------------------------------------------

df2_A11 <- df2 %>% filter(!is.na(CTxOficial))

m1 <- lm(Q.identidade_imp ~ CTxOficial, data = df2_A11)
m2 <- lm(Q.extr_lula      ~ CTxOficial, data = df2_A11)
m3 <- lm(Q.extr_bolso     ~ CTxOficial, data = df2_A11)
m4 <- lm(Q.human_lula     ~ CTxOficial, data = df2_A11)
m5 <- lm(Q.human_bolso    ~ CTxOficial, data = df2_A11)

stargazer(
  m1, m2, m3, m4, m5,
  type = "latex",
  title = "Table A11: Political Hostility (Official)",
  dep.var.labels = c(
    "Identity Salience",
    "Extremism Lulistas",
    "Extremism Bolsonaristas",
    "Humanization Lulistas",
    "Humanization Bolsonaristas"
  ),
  covariate.labels = c("CT X Official"),
  omit.stat = c("f", "ser", "adj.rsq")
)


# E.4 Robustness Checks --------------------------------------------------------

## Table A12 -------------------------------------------------------------------

summary(
  h2.elec.1.placebo <- lm(scale(Q.trust_eleições) ~ CTxPlacebo, df2)
)
summary(
  h2.elec.2.placebo <- lm(scale(Q.invasão) ~ CTxPlacebo, df2)
)
summary(
  h2.elec.3.placebo <- lm(scale(Q.aceitação1) ~ CTxPlacebo, df2)
)
summary(
  h2.elec.1.official <- lm(scale(Q.trust_eleições) ~ CTxOficial, df2)
)
summary(
  h2.elec.2.official <- lm(scale(Q.invasão) ~ CTxOficial, df2)
)
summary(
  h2.elec.3.official <- lm(scale(Q.aceitação1) ~ CTxOficial, df2)
)

h2.elec.1 <- lm(Q.trust_eleições ~ CTxOficial, data = df2)
h2.elec.2 <- lm(Q.aceitação1 ~ CTxOficial, data = df2)
h2.elec.3 <- lm(Q.invasão ~ CTxOficial, data = df2)

stargazer(
  h2.elec.1, h2.elec.2, h2.elec.3,
  type = "latex",
  dep.var.labels = c("Elec. Trust", "Elec. Accept", "8th January"),
  covariate.labels = c("Baseline: Official"),
  omit.stat = c("f", "ser", "adj.rsq"),
  digits = 3
)

## Table A13 -------------------------------------------------------------------

rm(list = ls())

# Load
att <- read.csv("attrition_df.csv", sep = ",", fileEncoding = "UTF-8")

# Attrition 
att$attrition <- ifelse(
  is.na(att$Q.extr_lula) | is.na(att$Q.extr_bolso) |
    is.na(att$Q.identidade_imp) | is.na(att$Q.human_lula) |
    is.na(att$Q.human_bolso) | is.na(att$Q.aceitação1) |
    is.na(att$Q.aceitação2),
  1, 0
)

# Groups
att$exp_arm <- dplyr::recode(
  as.character(att$Condition),
  "1" = "Right",
  "2" = "Left",
  "3" = "Freemasonry",
  "4" = "Placebo",
  "5" = "Official"
)

# Define baseline 
att$exp_arm <- factor(
  att$exp_arm,
  levels = c("Placebo", "Official", "Freemasonry", "Left", "Right")
)

# Model
attm <- lm(attrition ~ exp_arm, data = att)

# Table
stargazer(
  attm,
  type = "latex",
  title = "Table A15: Attrition Study 2",
  dep.var.labels = "Attrition",
  covariate.labels = c("Official", "Freemasonry", "Left", "Right"),
  omit = "Constant",
  omit.stat = c("adj.rsq", "f", "ser"),
  add.lines = list(
    c("Baseline", "Placebo")
  )
)

## Table A14 -------------------------------------------------------------------
## Table A15 -------------------------------------------------------------------

## Study 2 dataset
df2 <- read.csv('./Study 2_clean_pca outcomes.csv')

df2$exp_group <- factor(df2$exp_arm, levels = c('Placebo',
                                                'Official',
                                                'Freemasonry',
                                                'Left',
                                                'Right'))

df2$time <- ifelse(!is.na(df2$time_oficial_Page.Submit), df2$time_oficial_Page.Submit,
                   ifelse(!is.na(df2$time_placebo_Page.Submit), df2$time_placebo_Page.Submit,
                          ifelse(!is.na(df2$time_TCdireira_Page.Submit), df2$time_TCdireira_Page.Submit,
                                 ifelse(!is.na(df2$time_TCesquerda_Page.Submit), df2$time_TCesquerda_Page.Submit,
                                        ifelse(!is.na(df2$time_TCmaconaria_Page.Submit), df2$time_TCmaconaria_Page.Submit, NA)))))

df2 <- df2 %>%
  dplyr::mutate(complier = ifelse(exp_group == 'Placebo' & time > 112.75 |
                                    exp_group == 'Official' & time > 76.67 |
                                    exp_group == 'Freemasonry' & time > 126.9 |
                                    exp_group == 'Left' & time > 107.3 |
                                    exp_group == 'Right' & time > 121.25, 1, 0),
                attentive = ifelse(exp_group == 'Placebo' & time > 2*110 |
                                     exp_group == 'Official' & time > 2*73 |
                                     exp_group == 'Freemasonry' & time > 2*124 |
                                     exp_group == 'Left' & time > 2*104 |
                                     exp_group == 'Right' & time > 2*118, 0, 1))

# splitting time_video variable into bins of 3 seconds each
df2$time_video_bin <- as.numeric(cut2(df2$time,
                                      cuts = seq(min(df2$time),
                                                 max(df2$time), 3)))

# creating function to calculate statistical mode
getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

# creating continuous variable for compliance
# here we create a weighting variable based on the inverse of the absolute
# difference from the modal seconds bin to the actual seconds bin of each
# observation

# we sum the absolute difference by one not to divide by zero
df2 <- df2 %>%
  group_by(exp_arm) %>%
  dplyr::mutate(compliance_n = 1 / (abs(getmode(time_video_bin) - time_video_bin) + 1))


summary(
  h2.p.late <- lm(trust.pca ~ CTxPlacebo, df2 %>% filter(attentive == 1))
)

summary(
  h2.p.wls <- lm(trust.pca ~ CTxPlacebo, df2, weights = compliance_n)
)

summary(
  h2.o.late <- lm(trust.pca ~ CTxOficial, df2 %>% filter(attentive == 1))
)

summary(
  h2.o.wls <- lm(trust.pca ~ CTxOficial, df2, weights = compliance_n)
)

summary(
  h3.p.late <- lm(conflict.pca ~ CTxPlacebo, df2 %>% filter(complier == 1))
)

summary(
  h3.p.wls <- lm(conflict.pca ~ CTxPlacebo, df2, weights = compliance_n)
)

summary(
  h3.o.late <- lm(conflict.pca ~ CTxOficial, df2 %>% filter(complier == 1))
)

summary(
  h3.o.wls <- lm(conflict.pca ~ CTxOficial, df2, weights = compliance_n)
)


stargazer::stargazer(h2.p.late, h2.p.wls, h2.o.late, h2.o.wls, type = 'latex',
                     dep.var.labels = c('Institutional Trust'),
                     covariate.labels = c('CT'),
                     column.labels = c("(LATE)", "(WLS)", "(LATE)", "(WLS)"),
                     column.separate = c(1, 1, 1, 1),
                     omit = c('Constant'),
                     omit.stat = c('adj.rsq', 'f', 'ser'),
                     add.lines = list(c('Baseline', 'Placebo', 'Placebo', 'Official', 'Official')))

stargazer::stargazer(h3.p.late, h3.p.wls, h3.o.late, h3.o.wls, type = 'latex',
                     dep.var.labels = c('Partisan Hostility'),
                     covariate.labels = c('CT'),
                     column.labels = c("(LATE)", "(WLS)", "(LATE)", "(WLS)"),
                     column.separate = c(1, 1, 1, 1),
                     omit = c('Constant'),
                     omit.stat = c('adj.rsq', 'f', 'ser'),
                     add.lines = list(c('Baseline', 'Placebo', 'Placebo', 'Official', 'Official')))

################################################################################
## F Principal Component Analysis ----------------------------------------------
################################################################################

## ---------------------------- PCA models -----------------------------------##

library(tidyverse)
library(haven)
library(sjlabelled)
library(psych)
library(sjPlot)
library(Hmisc)
library(lavaan)
library(broom)
library(stargazer)

# Data --------------------------------------------------------------------

rm(list = ls())

## Study 1 data
df1 <- read.csv('./Study 1_clean.csv')

## Study 2 data
df2 <- read.csv('./Study 2_clean.csv')

# Overall support for democracy PCA ---------------------------------------

## Generates scree plot, so we know how many components are in the data
splot1 <- df1 %>%
  dplyr::select(acc_results, 
                acc_results_cand,
                pol_influence, 
                strongman,
                ban_protests, 
                ban_critics, 
                violence,
                agression_pt,
                agression_bs, 
                sup_dem) %>%
  fa.parallel(fm = 'ols')


### The screeplot suggests three components based on the eigenvalues

## Cronbach's alpha with all the items
df1 %>%
  dplyr::select(acc_results, 
                acc_results_cand, 
                pol_influence, 
                strongman,
                ban_protests, 
                ban_critics, 
                violence, 
                agression_pt,
                agression_bs, 
                sup_dem) %>%
  alpha(check.keys = T)### standardized alpha = 0.65

## Estimates principal component model for an unidimensional construct with all 
## the items of the democratic attitudes battery
pca1 <- df1 %>%
  dplyr::select(acc_results, 
                acc_results_cand, 
                pol_influence, 
                strongman,
                ban_protests, 
                ban_critics, 
                violence, 
                agression_pt,
                agression_bs, 
                sup_dem) %>%
  pca(nfactors = 1, rotate = 'varimax')

pca1

### We see that the variable of support for a "strong man" as the president does
### not load onto the unidimensional factor

## Reestimates PCA without strongman
pca2 <- df1 %>%
  dplyr::select(acc_results, 
                acc_results_cand, 
                pol_influence, 
                ban_protests, 
                ban_critics, 
                violence, 
                agression_pt,
                agression_bs, 
                sup_dem) %>%
  pca(nfactors = 1, rotate = 'varimax', scores = T)

pca2

## Scree plot
pca2_scree <- df1 %>%
  dplyr::select(acc_results, 
                acc_results_cand, 
                pol_influence, 
                ban_protests, 
                ban_critics, 
                violence, 
                agression_pt,
                agression_bs, 
                sup_dem) %>%
  pca(nfactors = 9, rotate = 'varimax', scores = T)

pca2_scree

### Now all the items load significantly onto the resulting factor

## Save eigenvalues and create variance table
eigenvalues <- pca2_scree$values
variance_table <- data.frame(
  Component = paste0("PC", 1:length(eigenvalues)),
  Eigenvalue = eigenvalues,
  Proportion_of_Variance = pca2_scree$Vaccounted["Proportion Var", ],
  Cumulative_Variance = pca2_scree$Vaccounted["Cumulative Var", ]
)

print(variance_table)


## Table A16 ------------------------------------------------------------------

## Variance table with stargazer
stargazer(
  variance_table,
  type = "latex", 
  summary = FALSE,
  digits = 3,
  title = "Variance Explained by Principal Components",
  rownames = FALSE
)

## Create object to plot scree plot
scree_data <- data.frame(
  Component = seq_along(eigenvalues),
  Eigenvalue = eigenvalues
)


## Figure A4 ------------------------------------------------------------------

## Scree plot 
screeplot <- ggplot(scree_data, aes(x = Component, y = Eigenvalue)) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2, color = 'grey60') +
  labs(x = "Component", y = "Eigenvalue") +
  theme_bw() +
  scale_x_continuous(breaks = seq(1, 9, 1)) +
  theme(panel.grid = element_blank())

ggsave('./screeplot.png', screeplot, 'png', width = 4, height = 3)

## Chronbach's alpha of the new scale
df1 %>%
  dplyr::select(acc_results, 
                acc_results_cand, 
                pol_influence, 
                ban_protests, 
                ban_critics, 
                violence, 
                agression_pt,
                agression_bs, 
                sup_dem) %>%
  alpha(check.keys = T)
### standardized alpha = 0.69, better than the former model.

## Checks the unidimensionality score
df1 %>%
  dplyr::select(acc_results, acc_results_cand, pol_influence,
                ban_protests, ban_critics, violence, agression_pt,
                agression_bs, sup_dem) %>% 
  unidim()

## Loadings table
loadings_table <- as.data.frame(unclass(pca2$loadings))
print(loadings_table)

## Save information for loadings plot
loading_df <- data.frame(
  Variable = rownames(loadings_table),
  Loading = loadings_table$PC1
)


loading_df$Variable[loading_df$Variable == 'acc_results'] <- 'Accept. Election Results'
loading_df$Variable[loading_df$Variable == 'acc_results_cand'] <- 'Candidates accept election results'
loading_df$Variable[loading_df$Variable == 'pol_influence'] <- 'Undue political influence'
loading_df$Variable[loading_df$Variable == 'ban_protests'] <- 'Ban protests'
loading_df$Variable[loading_df$Variable == 'ban_critics'] <- 'Ban critics'
loading_df$Variable[loading_df$Variable == 'violence'] <- 'Elec. Violence'
loading_df$Variable[loading_df$Variable == 'agression_pt'] <- 'Aggression Lulistas'
loading_df$Variable[loading_df$Variable == 'agression_bs'] <- 'Aggresion Bolsonaristas'
loading_df$Variable[loading_df$Variable == 'sup_dem'] <- 'Support for Democracy'


## Figure A5 ------------------------------------------------------------------

## Loadings plot
loadings <- ggplot(loading_df, aes(x = reorder(Variable, Loading), y = Loading)) +
  geom_bar(stat = "identity", fill = 'gray') +
  coord_flip() + # Horizontal bars for readability
  geom_hline(yintercept = 0.3, linetype = 2, color = 'gray60') +
  labs(x = "Variable", y = "Loading") +
  theme_bw() +
  theme(panel.grid = element_blank())

loadings

ggsave('./loadings.png', loadings, 'png', width = 3.5, height = 3.5)


# Institutional trust PCA -------------------------------------------------

## Scree plot with trust variables only
df2 %>%
  dplyr::select(Q.trust.parties : Q.trust.media,
                Q.trust_eleições, 
                Q.invasão, 
                Q.aceitação1) %>%
  fa.parallel()

## Despite the fact that we pre-registered hypotheses with an unidimensional
## index, the scree plot suggests two principal components. However, we still
## observe that most part of the variance is explained by the first component

## PCA with trust variables
pca3 <- df2 %>%
  select(Q.trust.parties : Q.trust.media, 
         Q.trust_eleições, 
         Q.aceitação1,
         Q.invasão) %>%
  pca(nfactors = 1, rotate = 'varimax')

pca3

### Support for the Jan 8th act does not load onto the institutional trust 
### construct. We drop it from the scale and perform the PCA again

pca4 <- df2 %>%
  select(Q.trust.parties : Q.trust.media, 
         Q.trust_eleições, 
         Q.aceitação1) %>%
  pca(nfactors = 1, rotate = 'varimax')

pca4

### Now all the items have good factor loadings

## Loadings table
loadings_table_trust <- as.data.frame(unclass(pca4$loadings))
print(loadings_table_trust)

## Save information for loadings plot
loading_df_trust <- data.frame(
  Variable = rownames(loadings_table_trust),
  Loading = loadings_table_trust$PC1
)

loading_df_trust$Variable[loading_df_trust$Variable == 'Q.trust.parties'] <- 'Trust Parties'
loading_df_trust$Variable[loading_df_trust$Variable == 'Q.trust.congress'] <- 'Trust Congress'
loading_df_trust$Variable[loading_df_trust$Variable == 'Q.trust.judiciary'] <- 'Trust Judiciary'
loading_df_trust$Variable[loading_df_trust$Variable == 'Q.trust.police'] <- 'Trust Police'
loading_df_trust$Variable[loading_df_trust$Variable == 'Q.trust.media'] <- 'Trust Media'
loading_df_trust$Variable[loading_df_trust$Variable == 'Q.trust_eleições'] <- 'Trust Elections'
loading_df_trust$Variable[loading_df_trust$Variable == 'Q.aceitação1'] <- 'Accept. Election Results'


## Figure A7 ------------------------------------------------------------------

## Loadings plot
loadings_trust <- ggplot(loading_df_trust, aes(x = reorder(Variable, Loading), y = Loading)) +
  geom_bar(stat = "identity", fill = 'gray') +
  coord_flip() + # Horizontal bars for readability
  geom_hline(yintercept = 0.3, linetype = 2, color = 'gray60') +
  labs(x = "Variable", y = "Loading") +
  theme_bw() +
  theme(panel.grid = element_blank())

loadings_trust

ggsave('./loadings_trust.png', loadings_trust, 'png', width = 4.5, height = 4)

## PCA for scree plot
pca4_scree <- df2 %>%
  select(Q.trust.parties : Q.trust.media, 
         Q.trust_eleições, 
         Q.aceitação1) %>%
  pca(nfactors = 7, rotate = 'varimax')

pca4_scree

## Save eigenvalues and create variance table
eigenvalues_trust <- pca4_scree$values
variance_table_trust <- data.frame(
  Component = paste0("PC", 1:length(eigenvalues_trust)),
  Eigenvalue = eigenvalues_trust,
  Proportion_of_Variance = pca4_scree$Vaccounted["Proportion Var", ],
  Cumulative_Variance = pca4_scree$Vaccounted["Cumulative Var", ]
) 

print(variance_table_trust)


## Table A17 ------------------------------------------------------------------

## Variance table with stargazer
stargazer(
  variance_table_trust,
  type = "latex", 
  summary = FALSE,
  digits = 3,
  title = "Variance Explained by Principal Components",
  rownames = FALSE
)

## Create object to plot scree plot
scree_data_trust <- data.frame(
  Component = seq_along(eigenvalues_trust),
  Eigenvalue = eigenvalues_trust
)


## Figure A6 -----------------------------------------------------------------

## Scree plot 
screeplot_trust <- ggplot(scree_data_trust, aes(x = Component, y = Eigenvalue)) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2, color = 'grey60') +
  labs(x = "Component", y = "Eigenvalue") +
  theme_bw() +
  scale_x_continuous(breaks = seq(1, 9, 1)) +
  theme(panel.grid = element_blank())

screeplot_trust

ggsave('./screeplot_trust.png', screeplot_trust, 'png', width = 4, height = 3)


# Partisan hostility PCA --------------------------------------------------

### Scree plot suggests one principal component

## Performs a model that extracts a single component
pca5 <- df2 %>%
  select(extr.outgroup, 
         dehum.outgroup, 
         Q.identidade_imp, 
         Q.aceitação2) %>%
  pca(nfactors = 1, rotate = 'varimax')

pca5

### All variables have good loadings

## Cronbach's alpha
df2 %>%
  select(extr.outgroup, 
         dehum.outgroup, 
         Q.identidade_imp, 
         Q.aceitação2) %>%
  alpha()

### Standardized alpha = 0.36

## Loadings table
loadings_table_hostility <- as.data.frame(unclass(pca5$loadings))
print(loadings_table_hostility)

## Save information for loadings plot
loading_df_hostility <- data.frame(
  Variable = rownames(loadings_table_hostility),
  Loading = loadings_table_hostility$PC1
)

loading_df_hostility$Variable[loading_df_hostility$Variable == 'extr.outgroup'] <- 'Outgroup extremism'
loading_df_hostility$Variable[loading_df_hostility$Variable == 'dehum.outgroup'] <- 'Outgroup dehum.'
loading_df_hostility$Variable[loading_df_hostility$Variable == 'Q.identidade_imp'] <- 'Id. importance'
loading_df_hostility$Variable[loading_df_hostility$Variable == 'Q.aceitação2'] <- 'Elections reliable'


## Figure A9 -----------------------------------------------------------------

## Loadings plot
loadings_hostility <- ggplot(loading_df_hostility, aes(x = reorder(Variable, Loading), y = Loading)) +
  geom_bar(stat = "identity", fill = 'gray') +
  coord_flip() + # Horizontal bars for readability
  geom_hline(yintercept = 0.3, linetype = 2, color = 'gray60') +
  labs(x = "Variable", y = "Loading") +
  theme_bw() +
  theme(panel.grid = element_blank())

loadings_hostility

ggsave('./loadings_hostility.png', loadings_hostility, 'png', width = 4, height = 4)

## PCA for scree plot
pca5_scree <- df2 %>%
  select(extr.outgroup, 
         dehum.outgroup, 
         Q.identidade_imp, 
         Q.aceitação2) %>%
  pca(nfactors = 4, rotate = 'varimax')

pca5_scree

## Save eigenvalues and create variance table
eigenvalues_hostility <- pca5_scree$values
variance_table_hostility <- data.frame(
  Component = paste0("PC", 1:length(eigenvalues_hostility)),
  Eigenvalue = eigenvalues_hostility,
  Proportion_of_Variance = pca5_scree$Vaccounted["Proportion Var", ],
  Cumulative_Variance = pca5_scree$Vaccounted["Cumulative Var", ]
) 

print(variance_table_hostility)

## Table A18 -----------------------------------------------------------------

## Variance table with stargazer
stargazer(
  variance_table_hostility,
  type = "latex", 
  summary = FALSE,
  digits = 3,
  title = "Variance Explained by Principal Components",
  rownames = FALSE
)

## Figure A8 ----------------------------------------------------------------

## Create object to plot scree plot
scree_data_hostility <- data.frame(
  Component = seq_along(eigenvalues_hostility),
  Eigenvalue = eigenvalues_hostility
)

## Scree plot 
screeplot_hostility <- ggplot(scree_data_hostility, aes(x = Component, y = Eigenvalue)) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2, color = 'grey60') +
  labs(x = "Component", y = "Eigenvalue") +
  theme_bw() +
  scale_x_continuous(breaks = seq(1, 9, 1)) +
  theme(panel.grid = element_blank())

screeplot_hostility

ggsave('./screeplot_hostility.png', screeplot_hostility, 'png', width = 4, height = 3)


# Electoral trust PCA -----------------------------------------------------

## Reverses support for Jan 8th (higher values = lower support)
df2$Q.invasão <- max(df2$Q.invasão, na.rm = T) - df2$Q.invasão

pca6 <- df2 %>%
  dplyr::select(Q.trust_eleições, Q.invasão, Q.aceitação1) %>%
  psych::fa(nfactors = 1, rotate = 'varimax')

pca6


# Factor loadings graphs ---------------------------------------------------

## OSDI

### Renames relevant variables
model2 <- df1 %>%
  select(acc_results, 
         acc_results_cand, 
         pol_influence, 
         ban_protests, 
         ban_critics, 
         violence, 
         agression_pt,
         agression_bs, 
         sup_dem) %>%
  mutate(`Accept. Results (Individual)` = acc_results,
         `Accept. Results (Candidate)` = acc_results_cand,
         `Use Political Influence` = pol_influence,
         `Ban Protests` = ban_protests,
         `Ban Gvt Critics` = ban_critics,
         `Jstfy Pol. Violence (General)` = violence,
         `Jstfy Violence (PTista)` = agression_pt,
         `Jstfy Violence (Bolsonarista)` = agression_bs,
         `Support for Democracy (abstract)` = sup_dem) %>%
  select(`Accept. Results (Individual)` : `Support for Democracy (abstract)`) %>%
  pca(nfactors = 1, rotate = 'varimax')


### Renames component
colnames(model2$loadings) <- c('OSDI')

### Plots graph
fa.diagram(model2, main = '', digits = 2, marg = c(0, 0, 0, 0))


## Institutional trust

### Renames relevant variables
model3 <- df2 %>%
  dplyr::select(Q.trust.parties:Q.trust.media,
                Q.trust_eleições,
                Q.aceitação1) %>%
  dplyr::rename(
    `Trust Parties` = Q.trust.parties,
    `Trust Congress` = Q.trust.congress,
    `Trust Judiciary` = Q.trust.judiciary,
    `Trust Police` = Q.trust.police,
    `Trust Media` = Q.trust.media,
    `Trust Elections` = Q.trust_eleições,
    `Accept. Results (Candidate)` = Q.aceitação1
  ) %>%
  psych::pca(nfactors = 1, rotate = "varimax")

### Renames component
colnames(model3$loadings) <- c('Institutional Trust')

### Plots graph
fa.diagram(model3, main = '', digits = 2, marg = c(0, 0, 0, 0))


## Partisan hostility

### Renames relevant variables
model4 <- df2 %>%
  dplyr::select(extr.outgroup, 
         dehum.outgroup, 
         Q.identidade_imp, 
         Q.aceitação2) %>%
  dplyr::rename(`Outgroup extremism` = extr.outgroup,
         `Outgroup dehumanization` = dehum.outgroup,
         `Identity importance` = Q.identidade_imp,
         `Partisan threat` = Q.aceitação2) %>%
  pca(nfactors = 1, rotate = 'varimax')


### Renames component
colnames(model4$loadings) <- c('Partisan hostility')

### Plots graph
fa.diagram(model4, main = '', digits = 2, marg = c(0, 0, 0, 0))


# Exports data ------------------------------------------------------------

## Extracts pca results
df1$dem_index_pca <- pca1[['scores']]
df2$trust.pca <- pca4[['scores']]
df2$conflict.pca <- pca5[['scores']]
df2$elections.pca <- pca6[['scores']]

## Saves datasets
write.csv(df1, './Study 1_clean_pca outcomes.csv')
write.csv(df2, './Study 2_clean_pca outcomes.csv')

################################################################################
## G Treatment effects on disaggregated measures of partisan hostility
################################################################################

rm(list = ls())

df <- read.csv('./study 2.csv')

## Table A19 --------------------------------------------------------------------
###### hypothesis 2: political conflict/animosity

# reverses variables so they indicate more conflictive attitudes
df$threat <- max(df$Q.aceitação2, na.rm = T) - df$Q.aceitação2
df$dehuman_b <- (max(df$Q.human_bolso, na.rm = T) - df$Q.human_bolso) / max(df$Q.human_bolso, na.rm = T)
df$demuman_l <- max(df$Q.human_lula, na.rm = T) - df$Q.human_lula / max(df$Q.human_bolso, na.rm = T)

# perceived threat by other political groups
summary(
  h2a <- lm(data = df,
            threat ~ CTxPlacebo)
)

# Dehumanization - Lulistas:
summary(
  h2b <- lm(data = df,
            demuman_l ~ CTxPlacebo)
)

# Dehumanization - bolsonaristas:
summary(
  h2c <- lm(data = df,
            dehuman_b ~ CTxPlacebo)
)

# Extremism perception difference:
summary(
  h2d <- lm(data = df,
            extr_diff ~ CTxPlacebo)
)

# identity strength
summary(
  h2e <- lm(data = df,
            Q.identidade_imp ~ CTxPlacebo)
)


# 8th of jan riot:
summary(
  h2f <- lm(data = df,
            Q.invasão ~ CTxPlacebo)
)

#### table with results
stargazer::stargazer(h2a, h2b, h2c, h2d, h2e, h2f, type = 'latex',
                     dep.var.labels = c('Threat', 'Dehuman. \nLulistas', 
                                        'Dehuman. \nBolsonaristas', 'Perceived \nExtremism',
                                        'Importance ID', '8th Jan.'),
                     covariate.labels = c('Baseline: Placebo'),
                     omit.stat = c('adj.rsq', 'f', 'ser'))

################################################################################
## H Estimations with PCA outcomes
################################################################################

## Table A20 -------------------------------------------------------------------

rm(list = ls())

# PCA Models

## Study 1 dataset
df1 <- read.csv('./Study 1_clean_pca outcomes.csv')

## Study 2 dataset
df2 <- read.csv('./Study 2_clean_pca outcomes.csv')

## Changes reference level of the experimental group variable
df1$exp_group <- factor(df1$exp_group, levels = c('Placebo',
                                                  'Oficial',
                                                  'Maçonaria',
                                                  'Esquerda',
                                                  'Direita'))

df2$exp_group <- factor(df2$exp_arm, levels = c('Placebo',
                                                'Official',
                                                'Freemasonry',
                                                'Left',
                                                'Right'))

df1$CTxOficial <- df1$XCTxOficial

### OSDI -- Baseline: placebo
summary(
  h1_placebo_pca <- lm(data = df1,
                       scale(dem_index_pca) ~ CTxPlacebo)
)

### OSDI -- Baseline: official theory
summary(
  h1_official_pca <- lm(data = df1,
                        scale(dem_index_pca) ~ CTxOficial)
)

### OSDI -- Experimental group factor
summary(
  h1_groups_pca <- lm(data = df1,
                      scale(dem_index_pca) ~ exp_group)
)

### Trust -- Baseline: placebo
summary(
  h2_placebo_pca <- lm(data = df2,
                       scale(trust.pca) ~ CTxPlacebo)
)

### Trust -- Baseline: official theory
summary(
  h2_official_pca <- lm(data = df2,
                        scale(trust.pca) ~ CTxOficial)
)

### Trust -- Experimental arms
summary(
  h2_groups_pca <-lm(scale(trust.pca) ~ exp_arm,
                     df2))

nobs(h2_placebo_pca);nobs(h2_official_pca)

### Electoral Trust -- Baseline: placebo
summary(
  h2_elec_placebo_pca <- lm(scale(elections.pca) ~ CTxPlacebo,
                            df2))

### Electoral Trust -- Baseline: official theory
summary(
  h2_elec_official_pca <- lm(scale(elections.pca) ~ CTxOficial,
                             df2))

### Hostility -- Baseline: placebo
summary(
  h3_placebo_pca <- lm(data = df2,
                       scale(conflict.pca) ~ CTxPlacebo)
)

### Hostility -- Baseline: official theory
summary(
  h3_official_pca <- lm(data = df2,
                        scale(conflict.pca) ~ CTxOficial)
)

### Generates table
m_list <- list(
  h1_placebo_pca, h1_official_pca,
  h2_placebo_pca, h2_official_pca,
  h2_elec_placebo_pca, h2_elec_official_pca,
  h3_placebo_pca, h3_official_pca
)

stargazer::stargazer(
  m_list,
  type = "latex",
  omit.stat = c("rsq", "adj.rsq", "f", "ser"),
  dep.var.labels = c("OSDI", "Institutional Trust", "Electoral Trust", "Partisan Hostility"),
  keep = c("CTxPlacebo", "CTxOficial"),
  covariate.labels = c("Placebo", "Official Theory"),
  omit = "Constant"
)


################################################################################
## I Models with covariate adjustment
################################################################################

## Table A21 -------------------------------------------------------------------

## Models with covariate adjustment
df2$income <- as.numeric(as.factor(df2$Q.renda))
df2$education <- as.numeric(as.factor(df2$Q.education))
df2$sex <- as.numeric(as.factor(df2$Q.gender))
df2$age <- as.numeric(as.factor(df2$Q.age))

df1$CTxOficial <- df1$XCTxOficial

### OSDI -- Baseline: placebo
summary(
  h1_placebo_covs <- lm(data = df1,
                        scale(dem_index) ~ CTxPlacebo + income + education + sex +
                          age)
)

### OSDI -- Baseline: official theory
summary(
  h1_official_covs <- lm(data = df1,
                         scale(dem_index) ~ CTxOficial + income + education + sex +
                           age)
)

## Trust -- Baseline: placebo
summary(
  h2_placebo_covs <- lm(data = df2,
                        scale(institutional_trust) ~ CTxPlacebo+ income +
                          education + sex + age)
)

### Trust -- Baseline: official theory
summary(
  h2_official_covs <- lm(data = df2,
                         scale(institutional_trust) ~ CTxOficial + income +
                           education + sex + age)
)


### Electoral Trust -- Baseline: placebo
summary(
  h2.elec_placebo_covs <- lm(scale(electoral_trust) ~ CTxPlacebo + income +
                               education + sex + age,
                             df2))

### Electoral Trust -- Baseline: official theory
summary(
  h2.elec_official_covs <- lm(scale(electoral_trust) ~ CTxOficial + income +
                                education + sex + age,
                              df2))

### Hostility -- Baseline: placebo
summary(
  h3_placebo_covs <- lm(data = df2,
                        scale(partisan_hostility) ~ CTxPlacebo + income +
                          education + sex + age)
)

### Hostility -- Baseline: official theory
summary(
  h3_official_covs <- lm(data = df2,
                         scale(partisan_hostility) ~ CTxOficial + income +
                           education + sex + age)
)

m_list_2 <- list(
  h1_placebo_covs, h1_official_covs,
  h2_placebo_covs, h2_official_covs,
  h2.elec_placebo_covs, h2.elec_official_covs,
  h3_placebo_covs, h3_official_covs
)

stargazer(m_list_2,
          type = 'latex',
          keep.stat = c('n'),
          covariate.labels = c('Placebo',
                               'Official Theory',
                               'Income',
                               'Education',
                               'Sex',
                               'Age'),
          dep.var.labels = c('OSDI',
                             'Institutional Trust',
                             'Electoral Trust',
                             'Partisan Hostility'))

################################################################################
## J Treatment effect of each experimental arm
################################################################################

rm(list = ls())

## Study 1 dataset
df1 <- read.csv('./Study 1_clean_pca outcomes.csv')

## Study 2 dataset
df2 <- read.csv('./Study 2_clean_pca outcomes.csv')

## Changes reference level of the experimental group variable
df1$exp_group <- factor(df1$exp_group, levels = c('Placebo',
                                                  'Oficial',
                                                  'Maçonaria',
                                                  'Esquerda',
                                                  'Direita'))

df2$exp_group <- factor(df2$exp_arm, levels = c('Placebo',
                                                'Official',
                                                'Freemasonry',
                                                'Left',
                                                'Right'))

## Creates data.frame with outcome variables
outcomes <- df1[, c('acc_results',
                    'acc_results_cand',
                    'strongman',
                    'pol_influence',
                    'ban_protests',
                    'ban_critics',
                    'violence',
                    'sup_dem',
                    'agression_bs',
                    'agression_pt')] 

## Changes column names
names(outcomes) <- c('Accept election results',
                     'Candidates should accept election results',
                     'Support for strongmen',
                     'Undue political influence',
                     'Ban protests',
                     'Ban government critics',
                     'Electoral violence',
                     'Support for democracy',
                     'Aggression bolsonaristas',
                     'Aggresion lulistas')


## Table A22 --------------------------------------------------------------------

levels(df1$exp_group) <- c('Placebo', 'Official', 'Freemasonry', 'Left', 'Right')

h1_arms_p <- lm(scale(dem_index) ~ exp_group, df1 %>% filter(exp_group != 'Official'))
h2_arms_p <- lm(scale(institutional_trust) ~ exp_group, df2 %>% filter(exp_group != 'Official'))
h3_arms_p <- lm(scale(partisan_hostility) ~ exp_group, df2 %>% filter(exp_group != 'Official'))

h1_arms_o <- lm(scale(dem_index) ~ exp_group, df1 %>% filter(exp_group != 'Placebo'))
h2_arms_o <- lm(scale(institutional_trust) ~ exp_group, df2 %>% filter(exp_group != 'Placebo'))
h3_arms_o <- lm(scale(partisan_hostility) ~ exp_group, df2 %>% filter(exp_group != 'Placebo'))

stargazer(h1_arms_p, h1_arms_o, h2_arms_p, h2_arms_o, h3_arms_p, h3_arms_o,
          type = 'latex',
          dep.var.labels = c('OSDI', 
                             'Institutional Trust', 
                             'Partisan Hostility'), 
          covariate.labels = c('Freemasonry', 'Left', 'Right'),
          add.lines = list(c('Baseline:', 'Placebo', 'Official', 'Placebo', 
                             'Official', 'Placebo', 'Official')),
          keep.stat = c('n'))


################################################################################
## K The Role of Belief in Conspiracy Theories ---------------------------------
################################################################################

## Table A23 -------------------------------------------------------------------

df1$belief <- ifelse(!is.na(df1$Q21_v1), df1$Q21_v1,
                     ifelse(!is.na(df1$Q21_v2), df1$Q21_v2,
                            ifelse(!is.na(df1$Q21_v3), df1$Q21_v3,
                                   ifelse(!is.na(df1$Q21_v4), df1$Q21_v4,
                                          ifelse(!is.na(df1$Q21_v5), df1$Q21_v5, NA)))))

table(df1$belief)

summary(
  belief1 <- lm(scale(dem_index) ~ CTxPlacebo*belief,
                df1)
)

# election results
summary(
  belief2 <- lm(scale(election_results) ~ CTxPlacebo*belief,
                df1)
)

# checks and balances
summary(
  belief3 <- lm(scale(checks_bal) ~ CTxPlacebo*belief,
                df1)
)

# civil liberties
summary(
  belief4 <- lm(scale(civil_lib) ~ CTxPlacebo*belief,
                df1)
)

# violence index
summary(
  belief5 <- lm(scale(violence_index) ~ CTxPlacebo*belief,
                df1)
)

# support for democracy
summary(
  belief6 <- lm(scale(sup_dem) ~ CTxPlacebo*belief,
                df1)
)

stargazer::stargazer(belief1, belief2, belief3, belief4, belief5, belief6,
                     dep.var.labels = c('OSDI', 'Elections', 'C and B',
                                        'Liberties', 'Violence', 'Supp. Dem.'),
                     covariate.labels = c('CT', 'Belief', 'CT X Belief'), 
                     omit = c('(Constant)'),
                     omit.stat = c('adj.rsq', 'f', 'ser'),
                     add.lines = list(c('Baseline', 
                                        'Placebo',
                                        'Placebo',
                                        'Placebo',
                                        'Placebo',
                                        'Placebo',
                                        'Placebo'))
)

## Table A24 ------------------------------------------------------------------
summary(
  belief7 <- lm(scale(dem_index) ~ XCTxOficial*belief,
                df1)
)

# election results
summary(
  belief8 <- lm(scale(election_results) ~ XCTxOficial*belief,
                df1)
)

# checks and balances
summary(
  belief9 <- lm(scale(checks_bal) ~ XCTxOficial*belief,
                df1)
)

# civil liberties
summary(
  belief10 <- lm(scale(civil_lib) ~ XCTxOficial*belief,
                 df1)
)

# violence index
summary(
  belief11 <- lm(scale(violence_index) ~ XCTxOficial*belief,
                 df1)
)

# support for democracy
summary(
  belief12 <- lm(scale(sup_dem) ~ XCTxOficial*belief,
                 df1)
)

stargazer::stargazer(belief7, belief8, belief9, belief10, belief11, belief12,
                     dep.var.labels = c('OSDI', 'Elections', 'C and B',
                                        'Liberties', 'Violence', 'Supp. Dem.'),
                     covariate.labels = c('CT', 'Belief', 'CT X Belief'), 
                     omit = c('(Constant)'),
                     omit.stat = c('adj.rsq', 'f', 'ser'),
                     add.lines = list(c('Baseline', 
                                        'Official',
                                        'Official',
                                        'Official',
                                        'Official',
                                        'Official',
                                        'Official')))


################################################################################
## L Exploratory Hypotheses ----------------------------------------------------
################################################################################

rm(list = ls())

## Study 1 dataset
df1 <- read.csv('./Study 1_clean_pca outcomes.csv')

## Study 2 dataset
df2 <- read.csv('./Study 2_clean_pca outcomes.csv')

## Changes reference level of the experimental group variable
df1$exp_group <- factor(df1$exp_group, levels = c('Placebo',
                                                  'Oficial',
                                                  'Maçonaria',
                                                  'Esquerda',
                                                  'Direita'))

df2$exp_group <- factor(df2$exp_arm, levels = c('Placebo',
                                                'Official',
                                                'Freemasonry',
                                                'Left',
                                                'Right'))

# Figure A10 ---------------------------------------------------------------

rm(list = ls())

## Study 1 dataset
df1 <- read.csv('./Study 1_clean_pca outcomes.csv')

## Study 2 dataset
df2 <- read.csv('./Study 2_clean_pca outcomes.csv')

table(df2$partisanCTxPlacebo)

summary(
  dem.pts.p <- lm(scale(dem_index) ~ partisanCTxPlacebo, 
                  df1 %>% filter(partisanCTxPlacebo != 'Non-political CT'))
  
)

summary(
  dem.nonpts.p <- lm(scale(dem_index) ~ partisanCTxPlacebo, 
                     df1 %>% filter(partisanCTxPlacebo != 'Political CT'))
  
)

summary(
  trust.pts.p <- lm(scale(trust.pca) ~ partisanCTxPlacebo2, df2)
)

summary(
  trust.nonpts.p <- lm(scale(trust.pca) ~ freemCTxPlacebo, df2)
)

summary(
  hostility.pts.p <- lm(scale(conflict.pca) ~ partisanCTxPlacebo2, df2)
)

summary(
  hostility.nonpts.p <- lm(scale(conflict.pca) ~ freemCTxPlacebo, df2)
)

summary(
  dem.pts.o <- lm(scale(dem_index) ~ partisanCTxOfficial, 
                  df1 %>% filter(partisanCTxOfficial != 'Non-political CT'))
)

summary(
  dem.nonpts.o <- lm(scale(dem_index) ~ partisanCTxOfficial, 
                     df1 %>% filter(partisanCTxOfficial != 'Political CT'))
)

summary(
  trust.pts.o <- lm(scale(trust.pca) ~ partisanCTxOfficial, 
                    df2 %>% filter(partisanCTxOfficial != 'Non-political CT'))
)

summary(
  trust.nonpts.o <- lm(scale(trust.pca) ~ freemCTxOfficial, 
                       df2)
)

summary(
  hostility.pts.o <- lm(scale(conflict.pca) ~ partisanCTxOfficial,
                        df2 %>% filter(partisanCTxOfficial != 'Non-political CT'))
)

summary(
  hostility.nonpts.o <- lm(scale(conflict.pca) ~ freemCTxOfficial,
                           df2)
)


# Placebo baseline — Partisan CT
summary(
  dem.pts.p <- lm(
    scale(dem_index) ~ partisanCTxPlacebo, 
    df1 %>% filter(partisanCTxPlacebo != 'Non-political CT')
  )
)

# Placebo baseline — Non-partisan CT  
summary(
  dem.nonpts.p <- lm(
    scale(dem_index) ~ relevel(factor(partisanCTxPlacebo), ref = "Placebo"),
    df1 %>% filter(partisanCTxPlacebo != 'Political CT')
  )
)
# Official baseline — Partisan CT
summary(
  dem.pts.o <- lm(
    scale(dem_index) ~ partisanCTxOfficial, 
    df1 %>% filter(partisanCTxOfficial != 'Non-political CT')
  )
)

# Official baseline — Non-partisan CT
summary(
  dem.nonpts.o <- lm(
    scale(dem_index) ~ relevel(
      factor(partisanCTxOfficial),
      ref = levels(factor(df1$partisanCTxOfficial))[grepl("Official", levels(factor(df1$partisanCTxOfficial)), ignore.case = TRUE)][1]
    ),
    df1 %>% filter(partisanCTxOfficial != 'Political CT')
  )
)
dem.nonpts.p.td <- tidy(dem.nonpts.p) %>%
  mutate(outcome = 'Support for Democracy', 
         Baseline = 'Placebo',
         Treatment = 'Non-partisan CT')

dem.nonpts.o.td <- tidy(dem.nonpts.o) %>%
  mutate(outcome = 'Support for Democracy', 
         Baseline = 'Official Theory',
         Treatment = 'Non-partisan CT')

dem.pts.p.td <- tidy(dem.pts.p) %>%
  mutate(outcome = 'Support for Democracy', 
         Baseline = 'Placebo',
         Treatment = 'Partisan CT')

trust.pts.p.td <- tidy(trust.pts.p) %>%
  mutate(outcome = 'Institutional Trust', 
         Baseline = 'Placebo',
         Treatment = 'Partisan CT')

hostility.pts.p.td <- tidy(hostility.pts.p) %>%
  mutate(outcome = 'Partisan Hostility', 
         Baseline = 'Placebo',
         Treatment = 'Partisan CT')

dem.pts.o.td <- tidy(dem.pts.o) %>%
  mutate(outcome = 'Support for Democracy', 
         Baseline = 'Official Theory',
         Treatment = 'Partisan CT')

trust.pts.o.td <- tidy(trust.pts.o) %>%
  mutate(outcome = 'Institutional Trust', 
         Baseline = 'Official Theory',
         Treatment = 'Partisan CT')

hostility.pts.o.td <- tidy(hostility.pts.o) %>%
  mutate(outcome = 'Partisan Hostility', 
         Baseline = 'Official Theory',
         Treatment = 'Partisan CT')

dem.nonpts.p.td <- tidy(dem.nonpts.p) %>%
  mutate(outcome = 'Support for Democracy', 
         Baseline = 'Placebo',
         Treatment = 'Non-partisan CT')

trust.nonpts.p.td <- tidy(trust.nonpts.p) %>%
  mutate(outcome = 'Institutional Trust', 
         Baseline = 'Placebo',
         Treatment = 'Non-partisan CT')

hostility.nonpts.p.td <- tidy(hostility.nonpts.p) %>%
  mutate(outcome = 'Partisan Hostility', 
         Baseline = 'Placebo',
         Treatment = 'Non-partisan CT')

dem.nonpts.o.td <- tidy(dem.nonpts.o) %>%
  mutate(outcome = 'Support for Democracy', 
         Baseline = 'Official Theory',
         Treatment = 'Non-partisan CT')

trust.nonpts.o.td <- tidy(trust.nonpts.o) %>%
  mutate(outcome = 'Institutional Trust', 
         Baseline = 'Official Theory',
         Treatment = 'Non-partisan CT')

hostility.nonpts.o.td <- tidy(hostility.nonpts.o) %>%
  mutate(outcome = 'Partisan Hostility', 
         Baseline = 'Official Theory',
         Treatment = 'Non-partisan CT')


models2 <- bind_rows(dem.pts.p.td,
                     dem.pts.o.td,
                     dem.nonpts.p.td,
                     dem.nonpts.o.td,
                     trust.pts.p.td,
                     trust.pts.o.td,
                     trust.nonpts.p.td,
                     trust.nonpts.o.td,
                     hostility.pts.p.td,
                     hostility.pts.o.td,
                     hostility.nonpts.p.td,
                     hostility.nonpts.o.td) %>% filter(term != '(Intercept)')

models2$Baseline <- factor(models2$Baseline,
                           levels = c('Placebo', 'Official Theory'))

models2$outcome <- factor(models2$outcome,
                          levels = c('Support for Democracy', 
                                     'Institutional Trust',
                                     'Partisan Hostility'))

models2$term <- models2$Baseline

results2.graph <- ggplot(models2 %>% filter(term != '(Intercept)'), 
                         aes(x = term, y = estimate, color = Treatment)) +
  geom_point(position = position_dodge(0.5), size = 1.75) +
  geom_errorbar(aes(ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error), 
                width = 0, position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = estimate - 1.65*std.error, ymax = estimate + 1.65*std.error), 
                width = 0, position = position_dodge(0.5), size = 1.25) +
  geom_hline(yintercept = 0, linetype = 2, color = 'grey60') +
  labs(x = 'Baseline', y = 'Standardized Coefficients') +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = 'bottom') +
  facet_wrap(~outcome)

results2.graph

ggsave('./study 2_pt_nonpt results.png', results2.graph, 'png', width = 6, height = 3.75)

################################################################################
## M Deviations from PAPs   ----------------------------------------------------
################################################################################

# Figure A11 ---------------------------------------------------------------

df2$lulista <- ifelse(df2$Q.lula == 5, 'Lulista', 'Non-Lulista')
df2$bolsonarista <- ifelse(df2$Q.bolsonaro == 5, 'Bolsonarista', 'Non-Bolsonarista')

df1$lulista <- ifelse(df1$affect_lula_c == 'Affection', 'Lulista', 'Non-Lulista')
df1$bolsonarista <- ifelse(df1$affect_bs_c == 'Affection', 'Bolsonarista', 'Non-Bolsonarista')

bind_rows(
  ## lulistas
  tidy(
    lm(scale(dem_index) ~ rightCTxPlacebo, df1 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'OSDI',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(dem_index) ~ rightCTxOfficial, df1 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'OSDI',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(trust.pca) ~ rightCTxPlacebo, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Institutional Trust',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(trust.pca) ~ rightCTxOfficial, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Institutional Trust',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(conflict.pca) ~ rightCTxPlacebo, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Partisan Hostility',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(conflict.pca) ~ rightCTxOfficial, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Partisan Hostility',
           Group = 'Lulistas'),
  
  ### bolsonarista
  
  tidy(
    lm(scale(dem_index) ~ rightCTxPlacebo, df1 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'OSDI',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(dem_index) ~ rightCTxOfficial, df1 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'OSDI',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(trust.pca) ~ rightCTxPlacebo, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Institutional Trust',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(trust.pca) ~ rightCTxOfficial, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Institutional Trust',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(conflict.pca) ~ rightCTxPlacebo, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Partisan Hostility',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(conflict.pca) ~ rightCTxOfficial, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Partisan Hostility',
           Group = 'Bolsonaristas')
  
) %>% 
  mutate(Outcome = factor(Outcome, levels = c('OSDI', 'Institutional Trust', 'Partisan Hostility')),
         Baseline = factor(Baseline, levels = c('Placebo', 'Official Theory'))) %>%
  filter(term != '(Intercept)') %>%
  ggplot(aes(x = Baseline, y = estimate, color = Group)) +
  geom_point(position = position_dodge(0.5), size = 1.75) +
  geom_errorbar(aes(ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error), 
                width = 0, position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = estimate - 1.65*std.error, ymax = estimate + 1.65*std.error), 
                width = 0, position = position_dodge(0.5), size = 1.25) +
  geom_hline(yintercept = 0, linetype = 2, color = 'grey60') +
  labs(x = 'Baseline', y = 'Standardized Coefficients', color = '', 
       title = 'Treatment: Right CT') +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = 'right') +
  facet_wrap(~Outcome) -> rct_lulistas_bolsonarista

###### left-ct and bolsonaristas-non-bolsonaristas -- lulistas-non-lulistas

bind_rows(
  ## bolsonaristas
  tidy(
    lm(scale(dem_index) ~ leftCTxPlacebo, df1 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'OSDI',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(dem_index) ~ leftCTxOfficial, df1 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'OSDI',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(trust.pca) ~ leftCTxPlacebo, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Institutional Trust',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(trust.pca) ~ leftCTxOfficial, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Institutional Trust',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(conflict.pca) ~ leftCTxPlacebo, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Partisan Hostility',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(conflict.pca) ~ leftCTxOfficial, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Partisan Hostility',
           Group = 'Bolsonaristas'),
  
  ### non-bolsonaristas
  
  tidy(
    lm(scale(dem_index) ~ leftCTxPlacebo, df1 %>% filter(bolsonarista == 'Non-Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'OSDI',
           Group = 'Non-Bolsonaristas'),
  
  tidy(
    lm(scale(dem_index) ~ leftCTxOfficial, df1 %>% filter(bolsonarista == 'Non-Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'OSDI',
           Group = 'Non-Bolsonaristas'),
  
  tidy(
    lm(scale(trust.pca) ~ leftCTxPlacebo, df2 %>% filter(bolsonarista == 'Non-Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Institutional Trust',
           Group = 'Non-Bolsonaristas'),
  
  tidy(
    lm(scale(trust.pca) ~ leftCTxOfficial, df2 %>% filter(bolsonarista == 'Non-Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Institutional Trust',
           Group = 'Non-Bolsonaristas'),
  
  tidy(
    lm(scale(conflict.pca) ~ leftCTxPlacebo, df2 %>% filter(bolsonarista == 'Non-Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Partisan Hostility',
           Group = 'Non-Bolsonaristas'),
  
  tidy(
    lm(scale(conflict.pca) ~ leftCTxOfficial, df2 %>% filter(bolsonarista == 'Non-Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Partisan Hostility',
           Group = 'Non-Bolsonaristas'),
  
  ### lulistas
  tidy(
    lm(scale(dem_index) ~ leftCTxPlacebo, df1 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'OSDI',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(dem_index) ~ leftCTxOfficial, df1 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'OSDI',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(trust.pca) ~ leftCTxPlacebo, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Institutional Trust',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(trust.pca) ~ leftCTxOfficial, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Institutional Trust',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(conflict.pca) ~ leftCTxPlacebo, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Partisan Hostility',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(conflict.pca) ~ leftCTxOfficial, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Partisan Hostility',
           Group = 'Lulistas')
  
) %>% 
  mutate(Outcome = factor(Outcome, levels = c('OSDI', 'Institutional Trust', 'Partisan Hostility')),
         Baseline = factor(Baseline, levels = c('Placebo', 'Official Theory'))) %>%
  filter(term != '(Intercept)' & Group != 'Non-Bolsonaristas') %>%
  ggplot(aes(x = Baseline, y = estimate, color = Group)) +
  geom_point(position = position_dodge(0.5), size = 1.75) +
  geom_errorbar(aes(ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error), 
                width = 0, position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = estimate - 1.65*std.error, ymax = estimate + 1.65*std.error), 
                width = 0, position = position_dodge(0.5), size = 1.25) +
  geom_hline(yintercept = 0, linetype = 2, color = 'grey60') +
  labs(x = 'Baseline', y = 'Standardized Coefficients', color = '', 
       title = 'Treatment: Left CT') +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = 'right') +
  facet_wrap(~Outcome) -> lct_bolsonarista_lulista

# freemasonry treatment

bind_rows(
  ## bolsonaristas
  tidy(
    lm(scale(dem_index) ~ freemCTxPlacebo, df1 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'OSDI',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(dem_index) ~ freemCTxOfficial, df1 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'OSDI',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(trust.pca) ~ freemCTxPlacebo, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Institutional Trust',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(trust.pca) ~ freemCTxOfficial, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Institutional Trust',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(conflict.pca) ~ freemCTxPlacebo, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Partisan Hostility',
           Group = 'Bolsonaristas'),
  
  tidy(
    lm(scale(conflict.pca) ~ freemCTxOfficial, df2 %>% filter(bolsonarista == 'Bolsonarista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Partisan Hostility',
           Group = 'Bolsonaristas'),
  
  ### lulistas
  tidy(
    lm(scale(dem_index) ~ freemCTxPlacebo, df1 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'OSDI',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(dem_index) ~ freemCTxOfficial, df1 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'OSDI',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(trust.pca) ~ freemCTxPlacebo, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Institutional Trust',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(trust.pca) ~ freemCTxOfficial, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Institutional Trust',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(conflict.pca) ~ freemCTxPlacebo, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Placebo',
           Outcome = 'Partisan Hostility',
           Group = 'Lulistas'),
  
  tidy(
    lm(scale(conflict.pca) ~ freemCTxOfficial, df2 %>% filter(lulista == 'Lulista'))
  ) %>% 
    mutate(Baseline = 'Official Theory',
           Outcome = 'Partisan Hostility',
           Group = 'Lulistas')
  
) %>% 
  mutate(Outcome = factor(Outcome, levels = c('OSDI', 'Institutional Trust', 'Partisan Hostility')),
         Baseline = factor(Baseline, levels = c('Placebo', 'Official Theory'))) %>%
  filter(term != '(Intercept)') %>%
  ggplot(aes(x = Baseline, y = estimate, color = Group)) +
  geom_point(position = position_dodge(0.5), size = 1.75) +
  geom_errorbar(aes(ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error), 
                width = 0, position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = estimate - 1.65*std.error, ymax = estimate + 1.65*std.error), 
                width = 0, position = position_dodge(0.5), size = 1.25) +
  geom_hline(yintercept = 0, linetype = 2, color = 'grey60') +
  labs(x = 'Baseline', y = 'Standardized Coefficients', color = '', 
       title = 'Treatment: Freemasonry CT') +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = 'right') +
  facet_wrap(~Outcome) -> fct_bolsonarista_lulista


bs.ls <- ggpubr::ggarrange(rct_lulistas_bolsonarista, lct_bolsonarista_lulista,
                           fct_bolsonarista_lulista,  ncol = 1,
                           common.legend = F, legend = 'bottom')

bs.ls

ggsave('./treatments_lulistas_bolsonaristas.png', bs.ls, 'png',
       width = 6, height = 10)

################################################################################
## N Full Regression Models ----------------------------------------------------
################################################################################

## Tables (A27 to A33) in the main script file. 
